function [Ahatmat, thetaipw, QTECVSelectIdx, ATECVSelectIdx, xi] = CVScoreCalType2(Y,X,A,taulist,xi,B)
    % The selection of basis set for the data with minimum cross-validation
    % value.
    % Inputs:
    %       Y:  the outcome variable Y_i in the data {Y_i, X_i, A_i}, for i =
    %           1,... 2*n    
    %       X:  the basic covariates X_i in the data {Y_i, X_i, A_i}, for i =
    %           1,... 2*n
    %       A:  the basic covariates X_i in the data {Y_i, X_i, A_i}, for i =
    %           1,... 2*n
    %       taulist: the sequence of quantile index used in the
    %       computations
    % Outputs:
    %       Ahatmat:        the computed score matrix, B*n*len(taulist)
    %       thetaipwmat:    the boostraped samples of thetaipw in ATE
    %       CVselectIdx:    the model indexes selected for treated and
    %       untreated group based on CV
    % 

    
    lentaulist      =   length(taulist);
    N               =   size(X,1);
    n               =   N/2;
    
    % Given the treatment status a = 0,1, compute the quantiles associated
    % with the tau in the taulist, and reconstruct it as a n*lentaulist
    % matrix for the following computation.
    Y1              =   Y(A == 1);
    Y0              =   Y(A == 0);
    Y1mat           =   repmat(Y1,1,lentaulist);
    Y0mat           =   repmat(Y0,1,lentaulist);
    
    Q1taulist       =   prctile(Y1,100*taulist);
    Q1taulistMat    =   repmat(Q1taulist,1,N/2)';
    Q0taulist       =   prctile(Y0,100*taulist);
    Q0taulistMat    =   repmat(Q0taulist,1,N/2)';
    
    % Here we construct the dependent variables for the following
    % coefficients computations
    % 1{Y(a,i) < q(a,tau)}, where q(a,tau) is the tau-th quantile for the
    % group with a = 0,1. 
    % Here is the B*lentaulist matrix Y1mat, Y0mat compare with the 
    % B*lentaulist quantile matrix Q1taulistMat, Q0taulistMat.
    A1CVmat         =   (Y1mat <= Q1taulistMat);
    A0CVmat         =   (Y0mat <= Q0taulistMat);       
    
    % Construct the basis for treated and untreated group
    X1              =   X(A == 1,:);
    BasisX1         =   BasisConstruct(X1);
    X0              =   X(A == 0,:);
    BasisX0         =   BasisConstruct(X0);

    CVmat0          =   nan(length(BasisX0),lentaulist);
    CVmat1          =   nan(length(BasisX1),lentaulist);
    for i = 1:length(BasisX0)
        % Given the basis matrix bX for treated group
        %            CVtheta1  = (bX'*bX)^(-1)*(bX'*A1mat)
        % here we only use the constructed dependent variable in treatment
        % group.
        % Regress dependent variables on each of the Basis above and compute the CV values
        % the computation is following the page 10 in 
        % "Nonparametric Sieve Regression: Least Squares, Averaging Least Squares, 
        %  and Cross-Validation".
        CVtheta0        =   (BasisX0{i}'*BasisX0{i})\(BasisX0{i}'*A0CVmat);
        CVem0           =   A0CVmat - BasisX0{i} * CVtheta0;
        CVhm0           =   diag(BasisX0{i}*((BasisX0{i}'*BasisX0{i})\BasisX0{i}'));
        CVmat0(i,:)     =   mean((CVem0.^2)./((1 - CVhm0).^2));

        CVtheta1        =   (BasisX1{i}'*BasisX1{i})\(BasisX1{i}'*A1CVmat);
        CVem1           =   A1CVmat - BasisX1{i} * CVtheta1;
        CVhm1           =   diag(BasisX1{i}*((BasisX1{i}'*BasisX1{i})\BasisX1{i}'));
        CVmat1(i,:)     =   mean((CVem1.^2)./((1 - CVhm1).^2));
    end
%     
    QTECVSelectIdx         =   nan(2,lentaulist);
    for tauIdx =  1:lentaulist
        QTECVSelectIdx(1,tauIdx)   = find(CVmat0(:,tauIdx) == min(CVmat0(:,tauIdx)),1);
        QTECVSelectIdx(2,tauIdx)   = find(CVmat1(:,tauIdx) == min(CVmat1(:,tauIdx)),1);
    end
   
    Ahatmat             =   nan(B,N,lentaulist);
    Basis               =   BasisConstruct(X);
    for tauIdx =  1:lentaulist        
        tmpbx0              =   Basis{QTECVSelectIdx(1,tauIdx)};
        tmpbx1              =   Basis{QTECVSelectIdx(2,tauIdx)};
        d0                  =   size(tmpbx0,2);
        d1                  =   size(tmpbx1,2);
        theta0              =   zeros(d0,B);
        theta1              =   zeros(d1,B);
        for b = 1:B
            Tmpxi           =   xi(b,:)';
            theta0(:,b)     =   ((repmat(Tmpxi,1,d0).*tmpbx0)'*tmpbx0/n)^(-1)*((repmat(Tmpxi,1,d0).*tmpbx0)'*A/n);
            theta1(:,b)     =   ((repmat(Tmpxi,1,d1).*tmpbx1)'*tmpbx1/n)^(-1)*((repmat(Tmpxi,1,d1).*tmpbx1)'*A/n);
            
            while (sum(isinf(theta1(:,b))) > 0) || (sum(isinf(theta0(:,b))) > 0) ...
                    || (sum(isnan(theta1(:,b))) > 0) || (sum(isnan(theta0(:,b))) > 0)
                
                Tmpxi           =   -log(1-rand(N,1));
                xi(b,:)         =   Tmpxi;
                theta0(:,b)     =   ((repmat(Tmpxi,1,d0).*tmpbx0)'*tmpbx0/n)^(-1)*((repmat(Tmpxi,1,d0).*tmpbx0)'*A/n);
                theta1(:,b)     =   ((repmat(Tmpxi,1,d1).*tmpbx1)'*tmpbx1/n)^(-1)*((repmat(Tmpxi,1,d1).*tmpbx1)'*A/n);
            end
        end   
        TmpAhat0            =   tmpbx0*theta0;
        Ahat                =   TmpAhat0;
        TmpAhat1            =   tmpbx1*theta1;
        Ahat(A == 1,:)      =   TmpAhat1(A == 1,:);
        Ahatmat(:,:,tauIdx) =   Ahat';
        
    end
    
    
    ATECVmat0              =   nan(length(BasisX0),1);
    ATECVmat1              =   nan(length(BasisX0),1);
    for i = 1:length(BasisX0)
        ATECVtheta0        =   (BasisX0{i}'*BasisX0{i})\(BasisX0{i}'*Y0);
        ATECVem0           =   Y0 - BasisX0{i} * ATECVtheta0;
        ATECVhm0           =   diag(BasisX0{i}*((BasisX0{i}'*BasisX0{i})\BasisX0{i}'));
        ATECVmat0(i)       =   mean((ATECVem0.^2)./((1 - ATECVhm0).^2));

        ATECVtheta1        =   (BasisX1{i}'*BasisX1{i})\(BasisX1{i}'*Y1);
        ATECVem1           =   Y1 - BasisX1{i} * ATECVtheta1;
        ATECVhm1           =   diag(BasisX1{i}*((BasisX1{i}'*BasisX1{i})\BasisX1{i}'));
        ATECVmat1(i,:)     =   mean((ATECVem1.^2)./((1 - ATECVhm1).^2));        
    end
    
    ATECVSelectIdx         =   [find(ATECVmat0 == min(ATECVmat0),1); ...
                                find(ATECVmat1 == min(ATECVmat1),1)];
    ATEbx0                 =   Basis{ATECVSelectIdx(1)};
    ATEbx1                 =   Basis{ATECVSelectIdx(2)};
    d0                     =   size(ATEbx0,2);
    d1                     =   size(ATEbx1,2);
    ATEtheta0              =   nan(d0,B);
    ATEtheta1              =   nan(d1,B);
    for b = 1:B
        Tmpxi             =   xi(b,:)';
        ATEtheta0(:,b)    =   ((repmat(Tmpxi,1,d0).*ATEbx0)'*ATEbx0/n)^(-1)*((repmat(Tmpxi,1,d0).*ATEbx0)'*A/n);
        ATEtheta1(:,b)    =   ((repmat(Tmpxi,1,d1).*ATEbx1)'*ATEbx1/n)^(-1)*((repmat(Tmpxi,1,d1).*ATEbx1)'*A/n);
    end
    Ahat0               =   ATEbx0*ATEtheta0;
    Ahat1               =   ATEbx1*ATEtheta1;
    Ahat                =   Ahat0;
    Ahat(A == 1)        =   Ahat1(A == 1);
    thetaipw            =   sum(xi'.*repmat(A.*Y,1,B)./Ahat)./(sum(xi'.*repmat(A,1,B)./Ahat)) -  ...
                            sum(xi'.*repmat((1-A).*Y,1,B)./(1-Ahat))./(sum(xi'.*repmat(1-A,1,B)./(1-Ahat)));
    
end
